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An inverse method to recover the star formation history and 
reddening properties of a galaxy from its spectrum. 
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Abstract. We develop a non-parametric inverse method to investigate the star formation rate, the metallicity 
evolution and the reddening properties of galaxies based on their spectral energy distributions (SEDs). This ap- 
proach allows us to clarify the level of information present in the data, depending on its signal-to-noise ratio (S/N). 
When low resolution SEDs are available in the ultraviolet, optical and near-IR wavelength ranges together, we 
conclude that it is possible to constrain the star formation rate and the effective dust optical depth simultaneously 
with a signal-to-noise ratio of 25. With excellent signal-to-noise ratios, the age-metallicity relation can also be 
constrained. 

We apply this method to the well-known nuclear starburst in the interacting galaxy NGC 7714. We focus on de- 
riving the SFR and the reddening law. We confirm that classical extinction models cannot provide an acceptable 
simultaneous fit of the SED and the lines. We also confirm that, with the adopted population synthesis models 
and in addition to the current starburst, an episode of enhanced star formation that started more than 200 Myr 
ago is required. As the time elapsed since the last interaction with NGC 7715, based on dynamical studies, is 
about 100 Myr, our result reinforces the suggestion that this interaction might not have been the most important 
event in the life of NGC 7714. 
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|> ! 1. Introduction 



The integrated spectra of galaxies contain information 
about the ages of their stellar populations, the metallicity 
of the stars and the effects of dust extinction. 

In order to infer the history of the star formation rate 
(SFR), one usually tries to match the spectral energy dis- 
tribution (SED) as closely as possible with model popula- 
tions computed with various scenarios for the SFR. The 
quality of the SED fit is assessed either qualitatively by 
visual inspection or quantitatively, e.g. by the x 2 ~test or 
a more general expression of the likelihood. Often, the 
SFR is represented by a number of discrete values or by 
a predefined analytic form. This approach, which could 
be called the direct or synthetic method, has been largely 
used in the field of population synthesis (see Heavens et 
al. 2000, Reichardt et al. 2001 for efficient recent imple- 
mentations). Since one usually keeps the number of free 
parameters as small as possible, the adopted functional 
forms for the SFR and age-metallicity relation (AMR) 
impose certain limitations on what type of model pop- 
ulations can be considered. Typically, these will be corn- 
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binations of instantaneous bursts and episodes of constant 
star formation. Some inverse methods (Craig and Brown 
1986, Tarantola and Valette 1982 a, b) deal with such a 
problem in an opposite way : one tries to determine the 
functional form of e.g. the SFR with as much freedom as 
possible, with a resolution in time that is dictated by the 
information contained in the data. Because of the latter 
property, these methods are called non-parametric. 

This work presents a non-parametric inverse method 
to estimate characteristics of galaxy evolution such as the 
SFR, the AMR and the intrinsic dust extinction. As for 
all such approaches, the method is based on a probabilis- 
tic formulation of inverse problems (in our case the for- 
malism of Tarantola & Valette 1982 a,b). We apply the 
method in the framework of evolutionary population syn- 
thesis. In other words, possible solutions are only sought 
among those compatible with our current understanding 
of star formation and evolution: the relative fraction of 
stars of various masses is not arbitrary but follows an 
initial mass function, and theoretical evolutionary tracks 
combined with stellar spectral libraries determine the pos- 
sible emission spectra of isochrone stellar populations. A 
probabilistic formulation for the alternative empirical pop- 
ulation synthesis has been developed recently in the para- 
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metric case by Cid Fernandes et al. (2001). As noted by 
these authors, the exploration of solutions to an inversion 
problem can be tackled as a minimization problem or with 
an adequate sampling algorithm for the space of param- 
eters. The second type of approach provides a complete 
description of the uncertainties on the estimated param- 
eters, but may become difficult to implement in practice 
when some of the unknowns are non-parametric functions 
of time, which can take an immense variety of shapes. 
Here a minimization procedure is adopted. In addition, 
specific tools are used in order to estimate the validity of 
the inverse procedure, like the a posteriori covariance and 
resolution, and the mean index. 

The nature of the available data determines many of 
the capabilities and limitations of inversion procedures. 
Sets of equivalent widths in the optical spectrum have 
been used because these measurements are relatively easy 
to acquire (Pelat 1997, Boisson et al. 2000). Cid Fernandes 
et al. (2001) combined equivalent widths and colours in or- 
der also to constrain dust extinction. However, they used a 
classical one parameter description of extinction and red- 
dening. The effects of dust are rarely that simple (Witt 
& Gordon 2000). The work we present here was partly 
motivated by previous studies of starburst galaxy spec- 
tra, which not only made it clear that average obscura- 
tion laws depend on the type of galaxy observed (Calzetti 
et al. 1994), but also showed that complex distributions 
of the stars and the dust in space can lead to significant 
local deviations from this average. Langon et al. (2001) 
studied the ~330 central parsecs of the nuclear starburst 
galaxy NGC 7714, and emphasized the effects of the dif- 
ferent optical depths of dust along various lines of sight 
within their small aperture. They demonstrated that in 
dusty objects it is necessary to combine SEDs and emis- 
sion lines from the optical, the ultraviolet and the near-IR 
spectral ranges if one aims at recovering the SFR over a 
wide range of ages and useful information on the wave- 
length dependent attenuation by dust. With those results 
in mind, we chose to apply the inversion method to data 
sets such as those of Langon et al. (2001): in this pa- 
per, the empirical constraints are combined low resolution 
SEDs in the three spectral ranges together with emission 
lines of Hn. Computation times were not prohibitive in 
the present case. 

The paper is organised as follows. Section 2 1| presents 
our model assumptions for the SEDs and reddening. 
Section [| introduces the inverse method. After applying 
this technique to simulated SEDs in Sect.||, we give new 
constraints on the SFR and the reddening law for the 
nuclear starburst of NGC 7714 in Sect.||. 



2. Modeling galaxy spectra. 

2.1. The basis B x (t) 

The intrinsic SED of a synthetic stellar population de- 
pends on the following model ingredients : 



— the SFR ip(t), which gives the total mass of stars of 
current age t, 

— the initial mass function (IMF), 

— the metallicity of the stars, Z(t), which reflects the 
metallicity of the interstellar medium from which the 
stars are born, 

— the results of stellar evolution calculations, 

— stellar atmosphere models, which give the flux at each 
wavelength for any star, 

— photoionization models, which add the nebular emis- 
sion to the SED. 

Without extinction, the simulated flux distribution Fx 
is written as follows: 



ti 



F x = / i>(t)B x {t,Z(t))dt 



(1) 



As ip{t) is always positive, it will be convenient to write 



ip(t) = ip exp(a(t)) 



(2) 



where ij) is a constant (arbitrary, but fixed). Bx(t, Z(t)) is 
the spectrum of a coeval stellar population of current age 
t. We considered ages between tt = 1 Myr and ti = 16 Gyr. 
This model basis implicitly assumes a certain initial mass 
function and rests on a particular set of stellar evolution 
and stellar atmosphere models. In the following, we will 
test the IMFs of Scalo (1998) and of Salpeter (1955), 
with a lower stellar mass limit of O.IMq and an upper 
mass limit of 120 M©. We construct basis spectra for a 
grid of ages and metallicities with the population synthe- 
sis code Pegase (Fioc & Rocca-Volmerange, 1997). They 
are based on the stellar evolutionary tracks of the Stellar 
Astrophysics group of the Astronomical Observatory of 
Padua (and extensions thereof), and on the semi-empirical 
spectral library of Lejeune et al. (1997, 1998). We refer to 
Fioc & Rocca-Volmerange (1997) for details. Time is sam- 
pled with logarithmic timesteps of ~0.1; the applications 
in the following sections will show this is fine enough to not 
limit the resolution of the estimated evolutionary prop- 
erties. The adopted grid in metallicity is Z=0.02, 0.016, 
0.012, 0.008. 

With extinction: 



Fx = / x/j(t)Bx(t,Z(t))f ext {X,t)dt 



(3) 



The wavelength dependence of f ext results from the as- 
sumed extinction model, as recalled below. 



2.2. Adopted extinction models. 

The extinction in galaxies depends on two factors : the 
spatial distribution of the dust and the properties of the 
dust grains, which we describe with an opacity coefficient 
kx ■ In Sects, [j.2.1 and 2.2.2 , we adopt simple models which 
assume the knowledge of the dust distribution and of kx ■ 
Opacity curves are known to vary from one galaxy to the 
other, or even from one line of sight to the other within 
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the Milky Way. For the foreground extinction due to our 
own Galaxy, we use the opacity curve of Seaton (1979) 
at optical wavelengths, and of Howarth (1983) in the UV 
and IR. This curve is an average over many lines of sight 
towards Milky Way stars. Due to the presence of graphite 
grains, it shows a characteristic feature around 2200 A. 
Because the dust in many galaxies has an extinction curve 
that lacks a 2200 A bump, it has been argued that the 
Small Magellanic Cloud (SMC) opacity is more convenient 
for the modeling of galaxies' intrinsic extinction (Gordon 
et al. 1997). For the intrinsic extinction of our synthetic 
galaxies, we adopt the SMC opacity curve of Prevot et al. 
(1984) as listed by Gordon et al. (1997), except where we 
explicitly allow it to vary. 

In reality, the interstellar medium of a galaxy is not ho- 
mogeneous, and is probably a mixture of different types of 
clouds (diffuse and compact) and different types of grains. 
Moreover, the light is scattered by the interstellar matter, 
so it is not possible to model the effective dust opacity in 
a simple way (e.g. Witt & Gordon 2000). This led Calzetti 
et al. (1994, 2000) to derive an obscuration curve empir- 
ically from the integrated spectra of a sample of diverse 
starburst galaxies (Sect. |2~2~3| ). In Sect. 2.2.4 we generalize 
this idea: it is suggested to use the inversion technique to 
recover both the wavelength dependence of the effective 
optical depth, and the time dependance of the SFR. 

2.2.1. Foreground dust screen. 

In this simplest model, one has: 



fext(\t) = e 



= exp(-0.921 fc A £(5 - V)(t)) (4) 



where E(B — V)(t) is proportional to the column density 
of dust. Stellar populations with different ages can be af- 
fected with different amounts of extinction, but all stars 
of a given age are affected identically. As already men- 
tioned, the values of k\ are those of Prevot et al. (1984). 
In this case, t\ is directly representative of the optical 
properties of dust grains on the line of sight towards stars. 
Representative plots of k\ can be found, for instance, in 
the review of Calzetti (2001). 

2.2.2. Dust clouds in front of the stellar population. 

In real galaxies, it is unlikely that all stars even of a 
given age see the same distribution of dust. Various mixed 
gas+stars models and clumpy gas distribution models 
have been considered in the literature, to allow for this 
natural complexity. The model considered here consists of 
dust clouds (or clumps) distributed between the stars and 
the observer. The number of clumps on the line of sight to 
a star of age t obeys Poisson statistics : n(t) describes the 
average number of clumps along the line of sight. With 
this model, a fraction exp(— n(t)) of all stars of age t are 
seen without any obscuration. 



As demonstrated in Appendix^, f ext has the following 
analytical expression for the clumpy model : 

f ext (X, t) = exp { -n(t) (1 - e- T ^) } = (5) 

In this equation, t AiC is the optical depth for a single cloud: 

T X . C = 0.921 k x E{B-V) c (6) 

E(B — V) c is the colour excess for one cloud. The result- 
ing extinction is described equivalently with an effective 
optical depth f\(t). 

Note that when t\ c « 1, f ext becomes indistinguish- 
able from a dust screen with total optical depth n(t) r\ iC , 
and when t A-c >> 1 the attenuation becomes independent 
of wavelength. In the former case, the attenuated spectra 
contain information on the product n(t) t AjC but the two 
factors cannot be separated; in the latter case, no signa- 
ture of extinction will appear in the spectrum (the dust 
would only be revealed by infrared emission) . 

2.2.3. Empirical extinction model of Calzetti et al. 

In the starburst extinction models of Calzetti et al. (1994, 
2000), / ea ;t(A,t) is an unknown which can be determined 
directly from observations assuming that all galaxies in 
the samples have similar star formation histories. 
Again, f ext is written: 



fext{\ t) 

where : 



f x (t)<xk x E(B-V)(t) 



(7) 



(8) 



E(B-V) is the colour excess to the stellar continuum. The 
empirical obscuration coefficient kx results from the wave- 
length dependence of the grain properties, as well as from 
the space distribution of the dust, averaged over a typi- 
cal starburst galaxy. The proportionality constant in the 
above relation is determined by Calzetti et al. (2000) from 
the energy budget of the galaxies, which requires far-IR 
observations. 

2.2.4. Effective optical depth. 

Inversion methods make it possible to explore whether the 
available spectrophotometric data contain enough infor- 
mation to recover both the time dependence of the SFR 
and the effective optical depth fx- Thus, we may consider 
the effective optical depth, fx, as an unknown. In this 
paper, we then restrict ourselves to the assumption of a 
constant fx in time. 

In fact, the opacity information lies principally in 
the emission lines. The differential optical depth between 
wavelengths Ai and A2 is given by: 



In 



(Ha 1 /Ha 2 )t? 



= TX 2 



(9) 



(H Al /H A2 ) 

where H A is the emission in the line at wavelength A. The 
subscripts R and correspond, respectively, to the red- 
dened and the intrinsic emission line ratios. In a similar 
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fashion, the slope of the spectrum yields complementary 
information on the opacity distribution. However, this ap- 
proach does not provide the ratio of absolute to differential 
extinction. If f\ is a solution, another possible solution is 
given by: 



fLtW = exp(-f A + A) 



(10) 



where A is an arbitrary constant. The solution for the SFR 
then becomes: 



i>'(t) =^(t)exp(-A) 



(11) 



So, in the framework described in this subsection, one 
only obtains information about relative fluctuations in the 
SFR. In practice, one could force the curve fx to vanish 
when A tends towards infinity or use the far-infrared 
thermal emission to determine the absolute extinction. 



3. The inverse method. 

The determination of scalar functions (such as the SFR 
and the extinction law) from observational data is an un- 
derdetermined problem, because the amount of observa- 
tional data is only finite and thus cannot provide the in- 
formation for every detail of these functions. Application 
of a straight inversion technique to Equation [| could be 
very sensitive to the noise in the data, and could well give 
mathematically correct but unphysical results (Craig and 
Brown, 1986). The problem must be regularized, which 
corresponds to a smoothing operation (Twomey, 1977; 
Tikhonov & Arsenine, 1976). 

3.1. A generalized least-squares approach. 

The inverse method we use comes from statistical tech- 
niques that have been applied in geophysical analyses 
(Tarantola and Valette, 1982 a,b; Tarantola & Nercessian 
1984; Nercessian et al. 1984). The method resembles 
Bayesian approaches in that a priori ideas about the un- 
knowns are used to regularize the inversion. 

The conditional probability density / p0 st(M|Z?) for the 
vector M of the unknown parameters, given the observed 
data D obeys: 



/ P oBt(M|£>) cx £(D|M)/ prior (M) 



(12) 



where C is the likelihood function and / pr i or stands for the 
a priori probability distribution for the model parameters. 
M will typically contain a sequence of values of the SFR at 
all ages, followed by the other adopted functional or dis- 
crete parameters (e.g. the metallicity at all ages, the opti- 
cal depth in the case of a dust screen, the effective optical 



depth at all wavelengths in the framework of Sect. 2.2.4). 
If we assume that both the a priori probability and the 
errors in the data are distributed as Gaussian functions, 
and that there are no correlations between data and model 



parameters, we can write 



f pos t(M\D) cx exp (- §(£> - g(M))* ■ C^ 1 ■ (D - g(M)) 



4(M — M )* ■ Cq 1 ■ (M — M ) 



(13) 



Mq is the mean of the prior (i.e. the a priori guess for 
the model parameters), and Cq is the variance-covariance 



matrix of the prior (see Sect. 3.2.1). The operator g as- 
sociates observable properties with model M: D — g(M) 
if the model parameters are correct and there is no noise 
in the data. In our case, Equation^ shows that g pro- 
vides the spectral energy distribution Fx through a set of 
equations of the form Fx — g\(ip, Z, ...). Cd is the data 
variance-covariance matrix, which describes the observa- 
tional uncertainties. The superscript * refers to the adjoint 
operator. 

The best estimate M minimizes the quantity : 



±{D-g(M))*-Cf-{D-g(M)) 

+\(M — Mq)* • Cq 1 ■ (M — Mq) 



(14) 



In the linear case, the minimum would be reached in one 
step (Tarantola and Valette, 1982 a, b): 



M = M + 
Co ■ G* ■ (Cd + G ■ Co ■ G*)~ 



(D-g(M )) 



(15) 



G is the matrix of partial derivatives of g. In our case 
however, g is not linear, and the minimum is reached it- 
eratively: 



M, 



[fe+i] 



M + 

Go ■ G* fc ] • (Cd -t- <jr[fc] • oo 
(D + G [k] ■ (M [k] - Mo) - g(M [k] )) 



Ik] ' °0 • Lr[fe]J 



(16) 



k counts the number of iterations and Gm is the matrix 



of partial Frechet derivatives at step k: 



G[k] 



= (°L) 

\9Mj [k] 



In subsequent equations we shall abbreviate 



S[k] Cd + • C ■ G 



[*] 



(17) 



(18) 



To control the algorithm convergence at step k, we test 
the stabilisation of x 2 , with : 



Xf fe] = (D-g(M [k] ))* -C d 



(D-g(M [k] )) 



(19) 



3.2. Choice of the model parameters and the prior Mq. 

3.2.1. The a priori variance-covariance operator Cq 

The model parameters may include single value parame- 
ters as well as functions (an example is described in detail 
in Appendix |B.2| ). For discrete parameters, Co is the ma- 
trix of the a priori variances and a priori covariances be- 
tween the parameters. Large values of the variances must 
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be chosen in the absence of any initial knowledge about 
the parameter values. Otherwise, the second term of ex- 
pression |l4] prevents the algorithm from moving away from 
a meaningless prior. Let us define a temporal varia ble u 
(u = log(i) or u = t for instance; see Appendix BT). For 
an unknown function such as a(u) (i.e. the star formation 
history of Equation |J), defined at each time point, Cp in - 
corporates a functional operator C a (see Equation B. 7). 
Given the limited number of data, we must assume some 
regularizing properties for a. 



C a (u,u') = o a (u)<7 a {u')CoT:(u,u') 



(20) 



The parameter a a is the prior standard deviation. It con- 
trols the acceptable amplitude of the deviations from the 
prior. Cor(ii, u') is the adopted autocorrelation function 
between two points of age u and u . We adopt: 



Cor(u, u') — exp 



(u - u 'f 

SO: 



(21) 



The parameter £ a is the prior correlation length. It 
defines the resolution (in age) with which one expects to 
recover the galaxy history functions. A large value of a a 
will allow for more contrast in the function a(u) than a 
small one. On the other hand, if too large a value of £ Q 
is chosen, the small (short timescale) details of the func- 
tion a will disappear and the parameter fluctuations will 
become too smoothed. The adequate choice of the two a 
priori parameters (£ Q and a a ) is conditioned by the sta- 
bility of the inversion process, which itself is determined 
by the extent of the available data and its signal-to-noise 
ratio (see also Sect. p.3.2] ). 

3.3. Quality of the inversion. 

3.3.1. The posterior variance-covariance matrix. 

Approximate internal errors on the estimation of the pa- 
rameters can be computed once the minimization algo- 
rithm has converged. This is done by considering the local 
behaviour of the posterior probability density distribution 
around the best solution (see also Reichardt et al. 2001, 
Moultaka & Pelat 2000). From a second-order expansion of 
the posterior density distribution in this neighbourhood, 
one obtains: 



(22) 



Again, Cm is a set of discrete values if the considered 
parameters have discrete values. In the functional case, 
Cm is a set of functions of two variables. For instance, 
if the star formation history is one of the unknowns, Cm 
incorporates terms of the form Ca/ iQ (u, u'). The diagonal 
terms such as Cm,o(u, u) measure how fast the exponent of 
Equation fL3| increases when a is moved away from its best 
estimate. Note that in general this a posteriori variance is 
lower than the a priori variance a 2 a of Equation The 
smaller the a posteriori variance is, compared to the a 
priori variance, the more we have increased our knowledge 
of the value of the parameter by using the available data. 



3.3.2. The resolving kernel and the mean index. 

Suppose that we knew the true model A/t ruc . Then the 
observed data would be (if the errors are negligible) : 



D = g{M tTU e) 



(23) 



After the linearisation of g near Mt rue , so that Gk = G = 
g, Equation[l5| shows how the deviation of the true model 
parameters from a nearby initial guess would be degraded 
into the derived one : 



M - M = C G*S- 1 G(M tTUC - M ) 
The operator defined by 
K := C G*S- 1 G 



(24) 



(25) 



describes this degradation of information and is called the 
resolving kernel. Let's again, for illustration, focus on the 
part of K relevant to the star formation history a(u). 
Equation^ writes: 



a(u) -a = K(u, u')(a tIUC (u') - a )du' 



(26) 



If the real star formation history deviated from a constant 
by a delta function at age us, the shape of the star forma- 
tion history obtained by the algorithm would be K(u, us). 
Because of the linearisation invoked, this response is only 
indicative. One of its major powers is to draw attention 
to potential degeneracies, that may appear in K(u, us) as 
very broad or multiple peaks. 

Another important and useful concept is a measure of 
the information present in the data. This is closely linked 
to the resolving kernel. Suppose that the parameter of 
interest is nearly constant within the width of the kernel 
K. Then Equation p4| gives 



(M-M )(«) = (M t , 



Mi 



) mean 



K(u,u')du' (27) 



The integral is called the mean index I(u) 

I(u) := ( K(u,u')du' (28) 



and has the following meaning. If I(u) has a very low value 
(<< 1) the model resulting from the minimization algo- 
rithm is expected to lie close to the prior, whatever this 
prior is. The quality of the M estimate thus is poor. But 
if I(u) ~ 1 the model obtained from the algorithm is close 
to the average true model. In other words, only I{u) ~ 1 
ensures that the available data contained significant in- 
formation on the estimated model parameter function. In 
practice, a compromise has to be found. A mean index 
close to 1 can be obtained by an increase in the smooth- 
ing length £ Q of the a priori variance-covariance operator. 
However, increasing ^ a reduces the resolution (in age). 
The adequate choice for £ Q is the smallest value that pro- 
duces a mean index close to 1. 
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Table 1. Priors for the inversions in Sect. | 



Dust clouds with individual E(B-V) 


=0.2 (Figs.^ andfy 


parameter 


prior 


prior a 


prior £ 


a(t) 





1 


0.3 - 0.5 [log(t)] 


n(t) 


1 


1 


0.3 - 0.5 jlog(t)] 


E(B-V) 


0.5 


0.25 


not applicable 


Z(f) 


0.014 


0.001 


0.5 [log(i)] 


Unknown effective optical depth (Fig. q). 


parameter 


prior 


prior a 


prior £ 


a(t) 





0.5 


0.7 [log(t)] 




1 


3(A<l^m) 


iooo A 






0.1 (A>l/zm) 





Note : a(t) is the natural logarithm of the SFR expressed in 
M yr- 1 . 



4. Determining SFR and reddening from mock 
spectra. 

We have generated mock spectra from a given SFR (a(t)), 
AMR (Z(t)), E(B-V) and average number of clouds (n(t)) 
in order to test and illustrate the inverse procedure. The 
spectra extend through parts of the UV [1180-1680 A], the 
optical [3690-6600 A and 7000-9000 A] and the near-IR 
[2.06-2.39 jttm], as appropriate for the subsequent analy- 
sis of NGC 7714 data (see Sect. ||). The spectral resolution 
matches that of the models, which are based on the library 
of stellar spectra of Lejeune et al. (1998; A/AA ~ 100). 
The IMF of Scalo (1998) is adopted. 

The noise added to the input spectra corresponds to 
a signal-to-noise ratio (S/N) equal respectively to 500 
(Fig. |l|) and 25 (Fig. [2]), and is assumed to be gaussian. 
Iterations were stopped when \ 2 was stable to one part 
in 10 4 (typically 7 to 9 iterations). The prior properties 
adopted for the figures are summarized in Table [|. The 
quality of the fits of the spectral energy distribution is 
evaluated by means of a reduced \ 2 > xl (the X 2 01 Eq. [l9| 
divided by the number of data points). The results are 
not affected much by the prior provided that the reduced 
X 2 and the mean index are close to 1. Details of the 
implementation of the inversion method are provided in 



Appendix B.2. 

In Fig. [I] the results obtained with S/N=500 are shown. 
The input parameters are drawn in full lines in the top 
panels. The input value of E(B-V) for individual clouds 
is 0.2. The dashed lines show the best estimates of these 
parameters resulting from the inversion. They led to a xt 
of 1.028. The dot-dashed lines show the standard devia- 
tions from these best estimates, computed from the di- 
agonal terms of the posterior variance- covariance matrix 
of Sect. |3.3.l[ For E(B-V), the resulting estimate is 0.194, 
and the posterior standard deviation is 0.03. 

The global features in the SFR are recovered well, 
whereas the smallest are not. The posterior resolution can 
be examined in the middle panels. They display the re- 
sponse of the inversion algorithm to delta-function inputs 
in SFR(i), n(t) and Z(t), located at ages log(i) = 6.8, 7.7, 
8.5 and 9.4. The computation is based on the resolving 



kernel K of Sect. |3~3~2| . It is seen that the resolution com- 
patible with the mock data is about 0.3-0.6 in log(age). 
The mean index for the SFR is uniformly close to 1 over 
the age range. Together with the fact that the posterior 
standard deviation is much smaller than the prior value, 
this tells us the SFR is well constrained. 

The average number of clumps and the AMR are not 
resolved as well as the SFR. The reddening of individ- 
ual clouds, E(B-V), has been recovered well. Note that its 
value is co mpatible with the recovery of n(t) (see end of 
Sect. 2.2.1 ). If E(B-V) had been found close to a degen- 
eracy with n(t) would have been suspected; if it had been 
very large both its own value and n(t) would have been 
very uncertain. 

Despite the high signal to noise ratio, the mean index 
associated with the AMR begins to weaken. Note that 
the smallest values of the mean index tend to appear at 
ages where the star formation rate is low: it is not sur- 
prising that little information on the mctallicity of the 
corresponding stars can be found in the data. 

Figure |^ shows how much information can be recov- 
ered when the available data has S/N=25. Only the global 
trends are recovered for the SFR and the average number 
of clumps. For the AMR, the resolution and the mean in- 
dex fall down after log(age)=7 (10 Myr) : the recovered 
AMR corresponds to the prior. Searching for the AMR 
should be avoided with this type of data. 

In Fig. ||, we test how well the method is able to re- 
construct the SFR and the effective optical depth f\ from 



simulated spectra with S/N=25 (see Appendix B. 3). Here, 
we do not consider the AMR as an unknown, because at 
this level of noise and spectral resolution, the information 



present in the data is too poor. As mentioned in Sect. 2.2.4 
we assume f\ is the same for populations of all ages. The 
three top panels show respectively the SFR, the resolution 
in age and the SFR mean index as in previous figures. Only 
a large smoothing length £ Q gives a decent mean index: the 
age resolution on the SFR is poor. 

The input and output effective optical depths are pre- 
sented in the bottom panel. The mean index for f\ takes 
values close to 1 in the UV, then drops progressively to 
poor values, of order 0.5 at 1 ^m and 0.1 around 2 ^m. The 
low index values at long wavelengths are due essentially to 
the degeneracy between the absolute v alues of the effective 
obscuration and of the SFR (see Sect. 2.2.4 ): with a small 
prior standard deviation, the prior sets the result (with 
larger prior standard deviations, the algorithm doesn't 
converge). In Fig.|| we have forced the near-IR match of 
f\, so input and output SFRs are comparable. As the UV 
emission can only come from young stars, the UV slope 
provides a strong constraint on that part of the reddening 
law. The emission of the young component at longer wave- 
lengths does not by itself match the data. The optical and 
near-IR fluxes directly provide a minimum contribution of 
older stars and the main constraints on their age distribu- 
tion. Freedom in the extinction law acts as a correction to 
improve the fit. A relatively long prior correlation length 
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Fig. 1. Inversion results from a spectrum with a signal-to-noise ratio of 500. The left column focuses on the the SFR, 
the middle column on the mean number of dust clouds on the line of sight, and the right column on the AMR. In 
the top panels, the solid lines describe the inputs used to construct the mock spectrum. The dashed lines are the best 
estimates from the noisy mock spectrum. The dot-dashed lines, and the middle and bottom panels allow us to asses 
the quality of the inversion. In the middle panels, the response functions K(u, us) are shown for ages us = log(i) = 
6.8 (solid), 7.7 (dashed), 8.5 (dot-dashed) and 9.4 (dotted). See text (Sect.|) for details. 



for f x (1000 A) avoids excessive freedom, that would lead 
to the fitting of the noise in the data. 

In view of the limited resolution, spectral coverage and 
signal-to-noise ratio assumed here, the results of Fig. H] are 
very satisfactory and the method is promising. 

5. An example : NGC 7714. 

NGC 7714 is a well studied interacting spiral galaxy which 
hosts a starburst in its nucleus (Arp 1966, Weedman ct al. 
1981, Calzetti 1997). The peaked ground based morphol- 
ogy of the starburst, the small inclination of the galaxy 
disk seen from our perspective, and the relatively low 
extinction suggested by the ultraviolet (UV) brightness 
(Markarian & Lipovetskij 1974), motivated detailed stud- 
ies of this apparently simple system. As a result of high 
resolution imaging and spectroscopy, the simplest models 
for the nucleus had to be successively excluded (Gonzalez- 
Delgado et al. 1995, 1999, and references therein). The 



central ~ 330 pc of the galaxy were found to contain a 
variety of stellar populations, ranging from very young 
(~ 5 Myr) to old (the population of the underlying spi- 
ral). Simple dust screen extinction models, that were ini- 
tially compatible with the scarce available data (Puxley 
& Brand, 1994), were found to be inappropriate (Lancon 
et al. 2001, hereafter LGLG01). 

The study of the central ~ 330 pc of the galaxy by 
LGLG01 was based on a "manual" exploration of a set 
of model star formation histories. The SFR was assumed 
to be a combination of a maximum of four components, 
each component being represented either by a standard 
spiral galaxy SFR, an instantaneous burst, an episode of 
constant star formation or an exponentially decreasing 
star formation episode. The exploration of this parame- 
ter space was guided by preliminary studies in individual 
wavelength ranges. Such an exploration method is tedious 
and cannot be exhaustive. A chance remains that the best 
fitting model may be missed. On the other hand, this ex- 
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ploration has provided a range of constraints that are con- 
sidered robust because they are common to most of the 
satisfactory model adjustments to the data. Among those 
were the following. 

(i) The nucleus of NGC 7714 has been forming stars 
off and on over the past several hundred millions of years 
at an average rate of the order of 1 M Q yr _1 , with a brief 
enhancement of a factor of a few about 5 Myr ago. 

(ii) The extinction even in the central 300 pc of 
NGC 7714 is inhomogeneous; for instance, most of the UV 
luminosity of this area is due to an obscuration-free line 
of sight towards a young cluster which does not coincide 
with the maximum Brackett7 (Br7) emission. 

(iii) Assuming a constant IMF in time, the recent star 
formation episodes have enhanced the stellar mass in the 
central area by at least 10%, and more likely by 25%. 

(iv) Most of the satisfactory models implied that 
the level of star formation had increased already several 
100 Myr ago. This timescale called for confirmation, as dy- 
namical modeling of the interacting system indicates that 
only about 100 Myr have elapsed since closest approach 
between NGC 7714 and NGC 7715 (Smith & Wallin 1992, 
Smith et al. 1997, Struck & Smith 2002). 



The existence of multiwavelength data and previous 
detailed investigations of plausible models for the central 
regions of NGC 7714 make this galaxy a target of choice 
for the application of our automatic inversion method. The 
available data are described in detail in LGLG01. The 
spectral coverage is as in Sect.^. We have degraded the 
resolution of the observations to match the wavelength 
sampling of the models. A correction for foreground ex- 
tinction due to the Milky Way was also applied (E(B- 
V)=0.08). The signal-to-noise ratio varies along the spec- 
trum with a typical value of 25. The gap in the optical part 
of the spectrum is due to the removal of data of poorer 
quality (due to telluric absorption). 



5.1. Estimation of the SFR with simple extinction 
models. 



In this section, we apply the three extinction models pre- 
sented previously in Sect. 2.2: a screen of SMC-type dust, 
clouds of SMC-type dust, and the empirical attenuation 
law of Calzetti et al. (2000). The priors are listed in Table|| 
and discussed in Sect. 5.2. 
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Fig. 3. Inversion results from a spectrum with a S/N = 
25, when the wavelength dependence of the effective opti- 
cal depth is also considered as an unknown. Metallicity is 
assumed to be known. The first three panels correspond to 
the first column of Figs.[l] and |^. The bottom panel shows 
the assumed optical depth (full line), the recovered one 
(dashed line) and the internal error (dot-dashed). 



The main results are plotted in Fig.y. The derived 
star formation histories and extinctions are in good global 
agreement with the previous study of LGLG01. A bimodal 
SFR is found, with two peaks respectively centered at 
log(age)=6.4 (2.5 Myr) and 8.5 (300 Myr). The average 
level of the SFR over the last few 10 8 yr is of the order 
of 1 M Q yr -1 . Stars younger than 10 9 yr contribute about 
20 % of the total stellar mass that has ever been produced. 
The most recent episode of activity however, despite its 
high star formation rate, only added a few percent to the 
stellar mass. 

The posterior standard deviation and mean index 
shows that the SFR is well constrained, except around 
log(age)=7.5. The stars at these ages contribute little to 
the light, compared to the youngest burst in the UV, or 
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different IMF slopes. 
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Scalo IMF 
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dust screen 
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dust clouds 


2.3 


2.0 


Calzetti 


2.1 


1.9 



to the bulk of the somewhat older intermediate age stars 
at optical wavelengths. The amount of obscuration is con- 
strained tightly for the youngest stars, and reasonably well 
for the oldest, but the mean index is poor at intermediate 
ages. 

In the case of the cloudy dust model, n(t) and the 
reddening per cloud are searched for. E(B — V) c = 0.17 is 
found for the individual clouds (with a posterior standard 
deviation of 0.03). n(t) varies between 1 and 2 at young 
and intermediate ages, and increases to about 6 for the 
old populations. n(t) and E(B — V) c are combined using 
Equation[| and || to provide the final colour excess plotted 
in Fig.[l| Note that this value is systematically smaller 
than the simple product n(t) E(B — V) c . 

The dependence of the colour excess on stellar age is 
consistent with the results of LGLG01. Relatively large 
values are found at old ages. This trend is necessary in 
order to reach the near-IR flux level without compro- 
mising the fit at optical wavelengths (Fig.|5j see however 
Sect. 

In the UV, the SMC opacity curve is the steepest of 
the reddening laws we have considered. With that law, 
the blue slope of the UV emission of NGC 7714 sets the 
strongest limitation on the amount of dust towards young 
stars. Consequently, a smaller SFR is sufficient to produce 
the observed level of UV flux. The shape of the opacity 
curves also explains the differences at old ages: the dif- 
ferent E(B-V) of Fig.|| produce old components of very 
similar shape and contributions in plots such as Fig.^. 

The reduced \ 2 values in Table|| show that the effect 
of the IMF on the quality of the best fits is small, the 



10 



J.-L. Vergely et al.: Inversion of galaxy spectra. 




og(age) 

Fig. 4. Properties obtained for the nucleus of NGC 7714, 
assuming a Salpeter IMF and solar metallicity (Z=0.02). 
The top panel shows the star formation histories. The mid- 
dle panel shows the differential extinction and the bot- 
tom panel the cumulative contributions of stars of various 
ages to the total stellar mass. Solid: Calzetti dust (the 
errors given by the posterior standard deviation are indi- 
cated with dotted lines in the top panel); dot-dashed: dust 
screen; dashed: dust clouds. 



Salpeter IMF being favoured. The clumpy dust model and 
the attenuation law of Calzetti et al. produce a slightly 
better agreement with the observed SED than a simple 
SMC dust screen. 

An important caveat for all the results of this section 
is that the Br7 emission line is underestimated by a fac- 
tor of two by the models, and this even though we do 
not apply stronger extinction to the emission lines than 
to the stellar continua (i.e. we deviate from the prescrip- 
tion of Calzetti et al. 2000 in this point). The missing Br7 
emission corresponds to a nebular emission excited by a 
young population that is not found by the inversion algo- 
rithm when the wavelength dependence of the extinction 
is one of the laws considered here. The young population is 
likely to be underestimated. We note that the same prob- 
lem was faced by LGLG01 when they sought to adjust 
only the continuum, and led them to invoke additional 
highly obscured young stars. A similar ad hoc fix would 
also be successful here. As a partial answer to the lack of 
Br7 photons, the age found by automated inversion for 
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Fig. 5. Adjustment of the NGC 7714 data obtained with 
the dust description of Calzetti. Thick line: data (gaus- 
sian fits to the emission lines have been subtracted be- 
fore resampling). Thin solid line: model SED (interpolated 
through regions with no data). Also shown are the re- 
spective contributions of stars in the following age ranges: 
< 10 7 yr (dotted), 10 7 - 10 9 yr (dashed), > 10 9 yr (dot- 
dashed) . 



the most recent starburst is younger than the 5 Myr of 
LGLG01. 

5.2. Discussion of the inversion assumptions 

It is important to analyse the effects of the adopted pri- 
ors. As expected, the numerical values of the priors given 
in Table ^ are of no consequence where the mean index is 
close to 1 , as long as they remain in a range in gross agree- 
ment with the flux level of the data, i.e. within reach of the 
iteration procedure given the adopted prior standard devi- 
ations. In practice, reasonable starting points are obtained 
after a small number of initially random trials which show 
into which direction the algorithm pulls. 

Increasing the prior standard deviations or reducing 
the correlation lengths significantly from the values in the 
table results either in loss of convergence, or in unaccept- 
ably small values of the mean indices, or in large poste- 
rior variances: the information in the data becomes insuf- 
ficient. 

As it is, the mean index of the extinction is low at 
intermediate ages. Nevertheless, the results are found to 
be quite stable against changes in the priors. In Fig.^, 
extreme results are plotted in the illustrative case of 
Calzetti-type dust. They were obtained by modifying the 
initial guesses of both the SFR and the extinction in order 
to pull towards higher or lower resulting SFRs. The shape 
of the SFR appears to be robust, although the mass ratio 
of old to relatively young stars (which is sensitive to small 
changes in the actual value of the early SFR) appears to 
be uncertain to within a factor of two. 
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smaller, and the intermediate age and young populations 
becomes more important in terms of total stellar mass. 
However, all the tests run here provide results within the 
global envelope of the results already discussed above or 
in the previous section. 

None of our attempts have a provided a solution 
in which the increase in the SFR responsible for the 
large intermediate age population occured less than about 
300 Myr ago. 




Fig. 6. Effects of the prior assumptions on the estimated 
parameters, in the case of Calzetti-type dust. The thick 
lines are the results of Fig.||. The thin solid line is the 
result of double iteration (see text). The two other lines 
are obtained with priors designed to pull the SFR towards 
low (dot-dashed) and high (dotted) values at intermediate 
ages (i.e. where the mean index of E(B-V) is low). The 
values of xt remain between 1.85 (double iteration) and 
1.94. 



Also plotted in F ig.|^ is the result obtained when the 
outputs of Sect. 5T are used as priors for a new itera- 
tion. This step provides a good test of the initial conver- 
gence. The changes are satisfactorily small. After success- 
fully performing the test in several cases, we have stopped 
applying it systematically. 

Finally, we have neglected until now the uncertain- 
ties due to possible aperture mismatch between the UV, 
optical and near-IR observations. LGLG01 estimate that 
errors in the relative flux levels of the three spectral seg- 
ments are below 15 %. Although it is possible to build this 
uncertainty directly into the inversion algorithm, we have 
more conservatively run inversions on a series of modi- 
fied data. Factors between 0.85 and 1.15 were applied to 
the UV and near-IR parts of the spectrum. The result- 
ing fits are poor (xt > 2-5) when the near-IR flux is in- 
creased above the reference value. The best fits (xl = 1-8) 
are obtained with a reduced near-IR flux. These solutions 
are appealing because they don't require higher extinc- 
tions for the old populations than for the intermediate age 
ones. The earliest star formation rates are correspondingly 



5.3. Discussion of star formation history 

Struck & Smith (2002) discuss possible causes for a star 
formation episode that would have started 300 Myr or 
longer ago, but favour none in particular. The strength 
and extent of the stellar ring of NGC 7714 and of the tidal 
tails of the system suggest that the galaxies are observed 
shortly after closest approach and well before their proba- 
ble merger (Toomre & Toomre 1972, Barnes & Hernquist 
1992, Gerber & Lamb 1994). What "shortly after" pre- 
cisely means remains dependent on the model. Struck & 
Smith (2002) caution that modest changes of orbital in- 
clination or distance of closest approach can significantly 
change the structure of the disk waves, and thus the cen- 
tral SFR in NGC 7714. At the same time, they admit a 
variety of imperfections in their final choices (for instance, 
they attribute the lack of a strong stellar ring in their hy- 
drodynamical simulation to a slightly overestimated im- 
pact parameter). Furthermore, the timescale of the for- 
mation and disruption of tidal features depends on impact 
speed (Dubinski et al. 1999). Although the choice of ini- 
tially parabolic galaxy orbits made by Struck & Smith is 
reasonable, the previous histories of the galaxies, or galaxy 
halo structures different from those adopted may have led 
to a slower encounter with longer timescales. 

Various dynamical or hydrodynamical simulations of 
mergers produce double star formation episodes over 
timescales of the order of 10 9 yr, but the second one usu- 
ally corresponds to the final merger, which is not appro- 
priate for our case study (Mihos & Hernquist 1994, Bekki 
1998, Gerritsen & Icke 1999). However, some models pre- 
dict that each of these SF episodes, especially the first one, 
will last for several 100 Myr (Bekki 1998). If these mod- 
els incorporated feedback mechanisms associated with star 
formation, such as local heating of the ISM that partially 
prevents star formation, the predicted SFR would be much 
more irregular (Gerritsen & Icke 1999) and might describe 
what happened in NGC 7714. 

From the simulations available in the literature, our 
conclusion is that the last interaction with NGC 7715 
might explain SF timescales of up to about 200 Myr, 
while longer SF timescales call for a previous event. A 
previous passage of the companion galaxy and related (or 
unrelated) bar instabilities, are options worth considering 
(Smith & Wallin 1992, Friedli & Benz 1995). 
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Fig. 7. Determination of the NGC 7714 star formation 
rate for two different metallicities, assuming a Salpeter 
IMF and the attenuation law of Calzetti et al. (2000). 
Solid line: Z=0.02; dashed line: Z=0.008. The layout is as 



in Fig.|. 



5.4. The influence of the assumed metallicity. 



Section 5.1 considers a solar metallicity for NGC 7714. 



Here we investigate the effect of the assumed metallicity 
on the estimated SFR. 

Fig. compares the SFRs obtained for Z=0.02 and 
Z=0.008. Due to the smaller opacities of metal poor stel- 
lar atmospheres, a stellar population of a given age is in- 
trinsically bluer at Z=0.008 than at Z=0.02. To reproduce 
the spectrum with a given effective extinction law, a shift 
of the star formation episodes to somewhat larger ages is 
therefore required at Z=0.008. Even the oldest stars are 
not red enough to explain the near-IR flux, and enhanced 
reddening towards the old populations results. The result- 
ing absorption is compensated with a higher early star 
formation rate, and the predominance of old stars in the 
mass budget is reinforced. 

The larger width of the recent star formation episode 
compensates for the lower value of the maximum SFR, 
and the total amount of UV light produced is similar to 
the Z=0.02 case. 



5.5. Using a new basis of model spectra. 

Spectrophotometric properties of stellar populations with 
complex star formation histories depend critically on the 
properties of the building blocks of evolutionary models, 
the so-called single stellar population models that define 
the basis B>. In addition to the basis of the previous sec- 
tions, we constructed a second set using the new grid of 
Mouhcine (2001) and Mouhcine & Lancon (2002). In this 
grid particular care was taken to model the near-infrared 
properties of intermediate age stellar populations. The 
major difference relative to the inputs of Pegase lies in 
the modeling of the Thermally Pulsing Asymptotic Giant 
Branch (TP-AGB). TP-AGB stars dominate the near- 
infrared light of intermediate-age (0.1-2. Gyr) stellar pop- 
ulations. As a result of the interplay between mass loss, the 
third dredge-up process and, most importantly, envelope 
burning, the new models predict that the TP-AGB stars 
with the longest lifetime are born from stars with main se- 
quence lifetimes of 0.8-1 Gyr rather than 0.1-0.2 Gyr (see 
also Girardi & Bertelli 1998, Marigo 1998). Consequently, 
the new integrated spectra evolve from bluer to redder 
colours between 0.1 and 1 Gyr rather than the opposite, 
in better agreement with constraints from the clusters of 
the Magellanic Clouds. One motivation for using the new 
models was to test whether they would allow for a more re- 
cent onset of the intermediate age star formation episode. 





og(oge) 

Fig. 8. Same as Fig.||, but with the stellar populations 
of Mouhcine & Langon (2002). Solid: Calzetti-type dust; 
dot-dashed: SMC dust screen; dashed: dust clouds. 
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Fig. U is the equivalent of Fig.|], for the new basis of 
single population spectra. The results are not significantly 
different. The intermediate age star formation episode 
starts early again. In the nucleus of NGC 7714, because 
of the large contribution of the old stars of the underlying 
spiral bulge, there is no need to invoke the reddest AGB- 
dominated populations that would differentiate between 
the two sets of basis spectra. 




s 
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5.6. Deriving the optical depth from the spectrum. 

The xt values of the previous models are of about 2, mean- 
ing that better models can be sought for. Moreover, the 
adjustment of the emission lines showed a specific prob- 
lem: Br7 is systematically underestimated by a factor of 
two. So, a possibility is that the extinction laws we have 
used do not correspond to reality. In this section, we derive 
a new extinction law, compatible with the data, emission 
lines included. More precisely, the model is written as fol- 
lows: 



Fx = (A) / </>o exp( tt (i)) B x {t, Z ) At 



with: 



/ext(A) = exp(-f A ) 



(29) 



(30) 



Here, we consider the effective optical depth f\ as an un- 
known but assume it is independent of age, and we adopt 
solar metallicity for all stellar populations. We apply the 
same extinction to the stellar continuum and to the nebu- 



lar emission. The basis spectra used are those of Sect. 5.1 
Appendix B.S provides computational details and Table [2| 
the priors. 

The SFR obtained here and shown in the top panel of 
Fig. ^| provides a good fit to the data (xl—lA; Fig.^Oj). 
The hydrogen lines are reproduced to within 10%. The 
mean index for the SFR is remarkably similar to the curves 
plotted in Figs.[I| or [| The mean index for f\ behaves as 
described and explained in Sect.^. The plotted attenua- 
tion law is derived to within an additive constant, which 
corresponds to multiplicative factor on the SFR. A mini- 
mum level of the SFR is set by the constraint of positive 
extinction at all wavelengths. If extinction were nil in the 
near-IR, the SFR would be reduced by only exp(— 0.1), 
i.e. about 10% as compared to the plotted curve. Much 
larger values of f\ and the SFR can be excluded when 
considering the far-IR emission of the galaxy, a constraint 
that we have not incorporated in the inversion procedure. 
LGLG01 verified that, with 4 times as many hot stars as 
directly derived from the UV flux with standard extinction 
laws, the nuclear ~ 330 pc would already produce over a 
third of the total far-IR light of the galaxy. This factor of 
4 is already present in the SFR of Fig.[| 

Interstingly, the wavelength dependence of the derived 
attenuation is close to the empirical law of Calzetti et 
al. (2000). The differences, however, are significant: in 
the UV, the estimated law is flatter (otherwise the blue 
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Fig. 9. Simultaneous determination of the star formation 
rate (top panel) and of the optical depth (middle panel, 
full line) for the central 300 pc of NGC 7714. We have su- 
perimposed the Calzetti law, rescaled to E(B — V) — 0.3 
(dashed line). The bottom panel shows the cumulative 
mass contribution of the stars. 

slope of the observed UV continuum could not be repro- 
duced with this amount of attenuation). The present re- 
sults mimic what imaging suggests is the real situation in 
the central 300pc of NGC 7714 (LGLG01): a significant 
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Fig. 10. Adjustment of the spectrum of the nucleus of 
NGC 7714, with the SFR and attenuation obtained in 
Sect. 5i. The line types are as in Fig.||. 
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fraction of the UV-emitting stars are heavily attenuated 
in the UV, and contribute mainly to Br7, while a smaller 
fraction of these hot stars is practically unobscured, ex- 
plaining the small apparent reddening at UV wavelengths. 

The SFR found here exceeds 1 M Q yr _1 only in the last 
20Myr. This timescale is consistent with the dynamical 
models of Struck & Smith (2002). Nevertheless, there is an 
enhancement in the SFR about 7 10 s yr ago (from about 
0.1 to about O.SMoyr" 1 ). 

Despite the appeal of this solution, we express some 
caution. First, the fit of Fig. [T^ shows the stars with ages 
below 10 7 yr as the strongest contributors around the 
Balmcr jump. The fit in this region is good because of 
a bump in the obscuration law at the appropriate place. 
There is no physical reason to expect such a bump. At 
higher spectral resolution, the shape of Balmer line wings 
would be able to confirm the need for a stronger contri- 
bution of intermediate age stars (LGLG01). Second, we 
have made the strong assumption here that a single em- 
pirical attenuation law applies to stellar populations of 
all ages. The solutions of LGLG01, based on the inte- 
grated spectrum and direct constraints from high reso- 
lution imaging, assigned populations of different ages dif- 
ferent apparent obscurations. In a way, we have traded 
freedom in the wavelength dependence of f\ against free- 
dom in its time dependence. Unfortunately, considering 
the quasi-infinite possible combinations of dust properties 
and dust distributions in a star burst, we doubt that any 
spatially integrated spectroscopic data will be able to fully 
break the degeneracies and make multiwavelength images 
redundant. 

6. Conclusions 

In this paper we have presented an inversion method de- 
signed to estimate the SFR, the AMR and the reddening 
properties of galaxies from their spectra. This method al- 
lows us not only to derive best values for these functional 
parameters, but also to estimate the amount of informa- 
tion actually present in the data for each of the unknowns, 
as well as the posterior resolution in time, depending on 
the S/N of the studied spectra. 

The main conclusions based on the inversion of simu- 
lated spectra are : 

— The SFR, the AMR and the reddening law are de- 
termined with a resolution of about 0.3-0.5 in log(age) 
for low spectral resolution UV + optical + near infrared 
spectra when noise in the data is negligible (S/N=500). 

— With S/N=25 and the same spectral coverage, the SFR 
is determined with a resolution of ^0.6 in log(age) 
and the information on the AMR is rather poor. 
Constraints on the reddening of the stellar continuum 
emission are tightest for the young stellar populations. 
The exact amount of recoverable information will de- 
pend on the actual star formation history of the ob- 
ject of study. Similar success has been reported by Cid 



Fernandes et al. (2001) in the framework of empirical 
population synthesis. 

— Alternatively, it is possible to constrain the SFR and 
the wavelength dependence of the effective optical 
depth simultaneously with a S/N=25, if it is assumed 
that this optical depth is independent of age (i.e. ap- 
plies to all stellar populations) and that the stellar 
metallicity is known. In that case, the resolution on 
the SFR is only about 1 in log(age). 

The inversion method applied to the spectrum of NGC 
7714 confirms that the reddening properties of this galaxy 
cannot be modeled with simple extinction models. In order 
to obtain satisfactory fits to both the stellar energy dis- 
tribution and the nebular emission lines, one has to allow 
the wavelength dependence of attenuation to deviate from 
"standard" laws (this work) or to allow for variations in 
the attenuation between different coexisting stellar popu- 
lations (LGLG01). 

The main results obtained for NGC 7714 can be com- 
pared to the previous results from non-automated studies 
(LGLG01): 

— The SFR found is indeed multimodal, with a young 
star formation peak a few Myr ago and a more ex- 
tended episode that formed stars at a typical rate of 
0.5-1 M yr _1 , thus providing an important interme- 
diate age population. 

— No solution was found in which the episode of en- 
hanced star formation that produced the bulk of the 
intermediate age stars started less than ~ 300 Myr 
ago. The discrepancy with the dynamical results of 
Smith & Wallin (1992), who estimate that ~ 100 Myr 
elapsed since closest encounter with NGC 7715, thus 
persists. New simulations tend to confirm the dynami- 
cal timescales (Struck & Smith 2002) , although a wider 
exploration of initial conditions may still allow it to be 
stretched a little. The importance of the event in terms 
of the stellar mass produced, especially when compared 
to the most recent starburst, tends to favour the oc- 
curence of an earlier perturbation, possibly unrelated 
with NGC 7715. 

— By solving simultaneously for the SFR and the effec- 
tive optical depth one strongly increases the estimated 
amplitude of the youngest peak in the SFR, while pro- 
ducing an attenuation law that is flattened in the UV 
part of the spectrum. The global shape of the attenua- 
tion is similar to the average law derived by Calzetti et 
al. (2000), but the differences in the details are signif- 
icant when the star formation rates are studied. The 
Br7 line emission is well matched, while it was under- 
estimated by a factor of 2 with "classical" attenuation 
laws. A drawback of this solution is that it produces 
a suspicious bump in the optical obscuration law in 
order to match the shape of the Balmer jump. 

— The solution just described mimics the effects that 
LGLG01 have shown to be due to a very inhomoge- 
neous distribution of the dust within the small region 
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observed, with one particular line of sight of very low 
extinction towards a cluster of hot stars. Even with 
UV, optical and near-IR spectroscopy combined, solv- 
ing for the effective optical depths requires strong as- 
sumptions. Here we have traded freedom in the wave- 
length dependence of optical depths against the free- 
dom to assign populations of different ages differing 
amounts of extinction. Nature might be even more 
complex, and require freedom in both when individ- 
ual galaxies are studied. 

The precious information of the spatial distribution 
of stars and dust within the region observed spectroscop- 
ically is not always accessible. The detailed study of ob- 
jects for which high resolution imaging is possible remains 
essential, and gives us insight into the fundamental uncer- 
tainties involved when only integrated light is available. 
Improvements will come naturally from a wider spectral 
coverage (if possible with overlaps between spectral ranges 
observed independently), and from the use of models with 
higher spectral resolution. 

Acknowledgements. We thank an anonymous referee for a care- 
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paper. 
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Appendix A: Dust clouds and Poisson law. 

When clouds are distributed homogeneously in space we 
can describe correctly the extinction as due to a Poisson 
distribution of clouds. The quantity n(t) represents the 
mean number of clouds encountered on the line of sight. 

We assume here that the clouds are of a typical size 
independent of age, so that the column density per cloud 
E(B — V) c is constant. In other words, we assume that 
the variation of the extinction in time depends only on 
the density of clouds. 

The probability to find i clouds in a given direction 
will follow a Poisson law : 

pU) = e -n(t)^fL ( A q) 
v. 

The proportion of the stellar population that is red- 
dened by i clouds is P(i). The absorption by i clouds 
is e -ikxE{B-v) c _ Then ^ the resulting flux ig (N a tta and 

Panagia 1984) : 

°° fti 

Fx = J2 mBx(t,Z(t))P(z)e- lk ^ B - v ^dt (A.2) 
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That yields : 



n(t) ST n< jy c -ik x E(B- 
i=0 ' 



with 

oc 



E 

i=0 



n W c -ik x E(B-V) c = e n(i)e- T A„ 



v ^dt 
(A.3) 

(A.4) 



The flux at the wavelength is modeled as follows : 

clog(t<) 



Appendix B: Estimation of the parameters. 

B.l. Change of variable: u = log(t). 

Massive stars evolve more rapidly than small mass stars, 
and a given absolute difference in initial masses has more 
effect on stellar evolution in the high mass regime than in 
the low mass one. As a consequence, the time interval over 
which spectrophotometric properties of isochrone stellar 
populations vary significantly increases quasi-linearly with 
age. The resolution in time expected for the SFR is thus 
almost constant in log(f) (see references in Hubeny et al. 
1999, and Lancon 2000). The following change of variable 
is natural: 



= log(t) 



(B.l) 



The fundamental equation (Equation^) should be written 
as follows: 

F x = / ^{u)B x {u,Z(u))f ext {\,u)h{u)du (B.2) 

Jlog(t f ) 

with h{u) = 10"ln(10). 

So, the prior variance-covariance operators defined be- 
low (C Q , Cz-, ...etc.) are expressed with the new variable 
it. For instance, C Q is written as : 

C a (u, v!) = a a (u)a a (u') exp I — 2 1 (B.3) 

so £ a is expressed in units of log(t). 

B.2. Determination of if), Z, n and E(B - V) 

In this section, we clarify the main steps of the inversion 
procedure, in the case one wishes to derive the time de- 
pendence of tp, Z and n, and the value of E(B — V) (the 
extinction per cloud on the line of sight). 

Since the star formation rate is always positive, chang- 
ing over from ip(u) to ce(u) in Equation^] yields: 



a(u) := ln(ijj(u)/tlJo) 



(B.4) 



where ipo is a constant. We shall assume that the uncer- 
tainties in the function a(u) follow a Gaussian law. This 
implies that the errors in ip(u) follow a log-normal dis- 
tribution. Similarly, the equations can be rewritten with 
logarithmic variables for all unknowns that must be forced 
to positive values. 



Fx, = / ip B Xi (u,Z{u)) 

J\og(t f ) 

exp {or(u) + n{u) (e^E^-V) _ A j h ^ du 

= g Xi (a(u),Z{u),n{u),E(B-V)) (B.5) 

This formulation requires use of the inverse method in the 
non-linear case (Equation fl6|). We drop the [k] index in 
order to simplify the notation. 
The vector M of unknown parameters is : 



M = 



( a{u) \ 

Z(u) 

n(u) 
V E(B - V) ) 



(B.6) 



The prior variance-covariance matrix shall be written : 
/C a \ 

r I Cz m 7 i 
c °- o o c w o (R7) 

V al [B _ v) ) 

Here C a (for instance) contains diagonal terms of the form 
C a (u,u) = er^(it) and non-diagonal terms of the form 
C a {u,u') (see Equation po|). 

The matrix of partial derivatives is (with s the number 
of wavelengths sampled): 



G 



\ 



dZ 
dg\ 2 
dZ 



dg\ s 
dZ 



9g\ ± dg Xl \ 



dg\ 



>><J,\ 



dE(B- 


v) 


dgx 2 




dE(B~ 


V) 


9gx s 





(B.8 



E(B — V) being a discrete parameter, the corre- 
sponding partial derivative is not a functional operator. 
dgx/dE(B - V) is written: 



dg 



x 



dE(B - V) 

%Bx(u,Z(u))n(u)kxe- k ^ B - V) 
. |a(«) + n(u) (e-kxE(B-v) _ ^ j h ^ du ^3 9^ 



-log(ti) 

/log( t/ ) 

exp • 



Conversely, a, Z and n are functions of the time vari- 
able. The corresponding partial derivatives are functional 
operators, ^p-, and have respectively the follow- 
ing kernels : 

K 1 (u)=^ B x (u,Z{u)) 

exp {a(u) + n{u) (e~^(B-v) _ X ) } h ^ ( B10 ) 



K 2 (u) = ip 



dB x (u,Z(u)) 
dZ 



exp {a(u) + n(u) [ e -^ E i B - v ) - l) } h{u) (B.ll) 
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K 3 (u) = [e 



-k x E{B-V) 



-1 i> B x (u,Z(u)) 



B.3. Determination oftp and f x . 



exp {a(u) + n(u) ( e - k \ E ( B ~ v ) — l) } h(u) (B.12) I n this section, the model and the observations are linked 

by the relation : 



These kernels Ki(u) act on a function f(u) as follows: 

"log(ii) 



<KiJ>-- 



f(u)K t (u)du 



(B.13) 



log(*/) 



exp(-f Ai ) 
g Xi (a{u),r x 



log(ii) 



ip B Xi (u, Zo)e a{u) h(u)dv(B.21) 



log(t/) 



where < , > is the scalar product in L 2 Hilbert space. 

The i-th component of the vector Vi = D + G ■ (M 
M ) - g{M) (in Equation^) is: 

dgx, 



where the metallicity is fixed. The unknown parameters 
are a(u) and f x with the respective variance-covariance 
operators C a and Cf • 



Vi = F\ 



dc 



-(a(u) - a ) 



M = 



a(u) 
f\ 



' gXt -(Z(u) - Z ) + ^^-{n{u) - no) 



dZ 



On 



dg Xl 



-(E(B-V)-E(B-V) ) 



c a 
c f 



dE(B - V) 
-g Xi (a(u),Z(u),n(u),E(B-V)) 



(B.14) 



Co = 
with : 

C f (A,A') = a T (\)a T (\')ex.p 



(A-A') 2 



(B.22) 
(B.23) 

(B.24) 



where ao,Zo, no and E(B — V)o are respectively the priors 
of a(u), Z(u),n(u) and E(B — V). The (i, j)-th component C a is given by Equation|2fj. To force positive values of f x , 



of (GC G*+C d ) is: 

c _ d 9^ n d 9*o , dgx z dg Xj dg Xz dg Xj 



the equations are rewritten with ln(rx) as the unknown. 
The matrix of partial derivatives of Equation is : 



da a da 



dE{B - V) 



dZ 



'E(B-V) 



dZ 
dg Xj 



dn 



dn 



dE(B - V) 



Si j <Ji<jj 



G = 



ogx 2 og\ 2 

da dr 



(B.15) 



da d 



(B.25) 



T I 



where b iti is the Kronecker symbol and a, the root mean ^ §m± respe ctively have the following kernels : 
square of the noise in the data F Xi (here we assume un- 
corrected noise in the data). 
Defining the vector : 



w = s~ l v, 



(B.16) 



Ki{u) = <tp eiq)(a{u)-T X )B x (u,Zo)h(u) (B.26) 
K 2 (\) = Tjjoexp(a(u)-T X )B x (u,Za)S(\-\i)h(u) (B.27) 

where S is the Dirac function. These kernels act on a func- 



tion as shown in Equation B.13 



the estimation number (fc + 1) for the parameters is com- 
puted from the previously estimated values by : 



a[k+i\{u) = a + ^2 W i 

i 

Z [k+1] (u) = Z Q + Y,W. l 



log(ti) 
\o g (t f ) 

log(ti) 
log(t f ) 

log(ti) 



C a (u, u')K\ (u')du' 



(B.17) 



C z {u,u')K 2 (u')du' 



The i-th component of the vector Vi = D + G ■ (M — 
M ) - g(M) is : 

Vi = F Xi + ——{a(u) - a ) + -^t-(t x - ta,o) 
da or 

-g Xi (a(u),f x ) (B.28) 

where ao and f Xi o are respectively the priors of a(u) and 
fx- 



(B.18) 



/•logUiJ 

i [k+1] (u) = n + V Wi / CW(u, u')K 3 (u')du' 

i ^log(i/) 



The matrix S (Equation B.15 ) is computed as follows: 

(B.29) 



c _ d 9Xi r 9g\i dgx t dg Xj 
da da ot x 8t x 



(B.19) 



where W is a vector defined by Equation B.16. The esti- 
mation of the parameters is given by : 



and 



E(B-V) [k+1] = E(B-V) 



U[k+i]{u) = a W i 



log(ti) 



l°g(*/) 



C a (u,u')Ki(u')du' (B.30) 



+ (J E(B-V)J2 m 



dg Xt 
dE(B - V) 



(B.20) 



f X[ k+i] = fx,o + Y,Wi Cf (A, X')K 2 (X')d\' (B.31) 
Jo 



